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ABSTRACT 

It is the purpose of this paper to develop a useful mathematical 
model of ASW aircraft availability. The increasing emphasis of systems 
studies dictates the use of accurate and representative models of the ASW 
systems. At present, many studies are using essentially the same models 
developed during World War II. This paper is an attempt to make use of 
advanced theory in a more powerful and flexible model and to make the use 
of the model practical and verifiable. 

The writer adapted the time homogeneous bivariate model as de- 
veloped by F. C. Collins. This is a discrete time Markov process with 
a stochastic matrix of transition probabilities wherein the maintenance 
process is modeled as a pulsed input multiple server queue. 

The model was programmed in FORTRAN 63 on the CDC l604and 
then modified to allow for variability, in the input parameters. Other 
modifications include an increase in the size of the model to accommodate 
a 1 6-aircraft squadron, the largest ASW squadron at present, and an 
explicit form solution to the maintenance queueing equations. 
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TABLE OF SYMBOLS 



Symbol 


Definition or Meaning 


a/ c 


aircraft 


\ 


mean repair rate of aircraft 


X A 


mean accident rate 


S 


the set of all possible outcomes; the probability 
description space 


E 


possible outcome (s) or event (s) 


(n, n + 1) 


a conditional probability that at time n + 1 the 
outcome or state is j given that at time n the 
state is i 


Xj (t) 


the number of a/c flying at time t, which did not 
fly the previous cycle 


X 2 (t > 


the number of a/c in the maintenance queue at 
time t 


A (t) 


the number of a/c desired on station at time t 


N (t) 


the total number of a/c of type considered at 
time t 


T 


the time interval from the launch to recovery at 
the start of the cycle 


Q(t o ) 


the probability distribution over all possible 
states at initial time t 

o 


P(t) 


the matrix of transition probabilities at time t 
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Symbol 



Definition or Meaning 



P (a, i) (P , 


the elements of the P matrix; the probability that 
X = (3 and X = j at the end of a cycle, given 

L C* 

that X = oi and X = i at the start of the cycle 


\h 


the probability given f ready a/c, g are launched, 
and h enter maintenance 


P Y 


probability of entering maintenance just before, 
during, or immediately after launch 


P 


the probability of equipment failure during flight 
requiring maintenance when recovered by the 
carrier 


n (m) 

a 


the probability that of a a/c flying m will enter 
maintenance upon recovery 


D 


the number of independent identical maintenance 
repair stations or "spots 11 


p y |l) 


the probability that i - j a/c are repaired in 
time interval t 
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1. INTRODUCTION 



The threat to freedom of the seas posed by the vast Soviet sub- 
marine fleet is perhaps the most thorny problem facing the U. S. Navy 
today. Two world wars have produced Pyrrhic victories over limited 
submarine fleets. During the Second World War operations analysis 
was born into the Navy to aid in the defeat of the German submarine. 

The classic antisubmarine warfare (ASW) analyses and models developed 
by Morse [2] and Koopmans [3] are still being used today, over two 
decades later, in most of the ASW study efforts for the Navy. 

These early ASW analyses assumed a given level of search effort 
available and directly evaluated the probability that an ASW subsystem 
could detect and/or kill a submarine. This assumption is not only 
logical to make the problem tractable, but also practical since no 
immediate changes in ASW force levels could be expected. Moreover, 
the studies were conducted during the war, not before it started. It is 
the purpose of this paper to present a probabilistic model to describe 
the available effort. Such a model can be used to sharpen the estimates 
of the effectiveness of an ASW subsystem and to study the charac- 
teristics of the associated support system. 

Naturally, the current study plays an important but limited role in 
the overall problem of designing an entire ASW system. The diffi- 
culties involved in such a specification are legion. First and foremost 
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is the quantification of the ASW mission in denying the enemy the 
effective use of his submarines. Currently, the probability of detecting 
and/or killing submarines is used as the measure of effectiveness of 
the mission, and it appears that a more encompassing one has not been 
developed. Second, the specification of an ASW force level to counter 
a given threat has many inherent subjective elements. These are due to 
the existing historical bias in predicting the conduct of a future ASW war 
with an enemy, particularly one who has never before used a large 
submarine force in its military operations. The reader can imagine 
why merely defining terms such as "threat 11 and "effective counter" 
becomes quite difficult. 

Thus, there is a need to investigate the levels of search effort 
specified. This may require acceptable models to measure the availa- 
bility of effort, its effectiveness, and determine the logistic support 
required for any level of available effort. Specifically, the ASW sub- 
system to be modeled is the carrier -based aircraft, although the model 
is adaptable to other systems. 

The method of investigating the demand for ASW carrier a/c will 
assume that the desired number of a/c on station is known as an input 
parameter. The support required to achieve this measure of available 
effort depends upon maintenance space, manpower, and supply. 
Generally, we shall consider how an ASW carrier supports this number 
of a / c on station with the present or proposed number of a/c embarked 
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on the carrier. The parametric input can be subjected to sensitivity 
analyses. 

The operational commander of the ASW force launches the desired 
number of a/c on station to screen, search, or actively prosecute a 
submarine contact. Each a/c is relieved on station. Each such relief 
requires the launching of another a/c prior to the recovery of the initial 
a/c. The returning a/c must receive varying degrees of maintenance 
and requires refueling and rearming. This cycle continues until the 
mission is completed. Loss of a/c due to accidents, insufficient supply, 
and lack of repair capability cause deviations in this procedure. Naval 
operations involve the interaction of many quantities which are random 
in nature. Not all can be considered in a tractable mathematical model. 
Some quantities which are important are omitted. One example is the 
length of each cycle time, which is assumed to be a constant value. 
Including variables of this nature incurs unnecessary mathematical 
complication. It is hoped that adequacy of the model can be measured 
by using fleet data available from the Fleet ASW Data Analysis Program 
(FADAP). 

Collins [5] describes a bivariate Markov model for airborne early 
warning (AEW) and combat air patrol (CAP) jet a/c operating in an 
attack carrier force. This model is used to evaluate the probability of 
maintaining a fixed requirement of a/c on station as a measure of ef- 
fectiveness of the system. It has subsequently been used in a larger 
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attack force study for the Navy. The model computes the probabilities 
of the number of a/c on station and in or awaiting maintenance at any 
given launch period. The comparable ASW problem differs in the following 
aspects: 

1. Type, range, and speed of a/c; 

Z. The variable number of a'/c required for mission; 

3. Attrition due to accidents and supply failures; 

4. The greater number of ASW a/c. 

It was decided to use the Collins' model with appropriate modification. 

For immediate reference, the mathematical content of the model will 
be repeated herein. 

In order to incorporate these modifications, it was necessary to 
spend some time reprogramming on the CDC 1604 digital computer in 
FORTRAN 63, the CDC version of the IBM FORTRAN IV. The original 
program was not readily available and was written in an early assembler 
language. Moreover, the numerical analysis was not sufficiently sharp 
to handle the larger input values. Also, double precision (two computer 
words instead of one) arithmetic was required in one subroutine for an 
accurate explicit solution to the maintenance queueing equations (see 
Appendix I). This effected a 50% decrease in the computer time 
required for developing a matrix of transition probabilities. 

Following this introduction, section Z contains a brief description of 
the operational problems involved and the assumptions made. A brief 
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description of Markov chains and the mathematical model are presented 
in section 3. The details for computing the matrix of transition proba- 

4 

bilities are given in section 4. General employment of the model 
follows. The appendices include the solution mentioned on the pre- 
ceding page, a logical flow diagram of the program, a copy of the 
program, and some sample results. 
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2. ASSUMPTIONS 



The real-world employment of carrier a/c is cyclic in nature, and 
the present state of any given a/c (i. e. , flying, in or awaiting mainte- 
nance) depends largely on what the previous state was. This fact 
suggests that a Markovian assumption can logically be made for the 
a/c transition probabilities . In the search phase, a/c may or may not 
relieve on station; but, in any part of the contact investigation phase, 
relief on station will be made. To insure full screening and mission 
coverage, a/c will relieve on station. 

The question of resupply during an operation depends primarily on 
the availability of carrier on-board delivery (COD). This depends on 
the geographical location and the mission (convoy protection, strike- 
force protection, hunter -killer operation, etc.). In practice, resupply 
is not anticipated within a week's period, and around-the-clock oper- 
ations have continued for two weeks without resupply. 

Standard maintenance procedures aboard carriers preclude major 
maintenance on the flight deck. It will be assumed that sufficient notice 
is given so that all major 120 -hour checks will be completed prior to 
the operation. This assumption can be modified with an appropriate 
adjustment in the mean repair rate. The concept of maintenance crews 
assigned to hangar deck areas ("spots"), as developed by Collins [3], 
will be used. Each crew will be capable of all types of maintenance 
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and will operate independently at the identical mean repair rate 
The number of spots is determined by the average number of such 
crews available to work continuously around the clock on a watch basis. 

The state of each a/c is assumed to be statistically independent of 
that of others, and the launching and landing transition probabilities 
will be developed on the basis of independent Bernoulli trials. The 
parameters can be determined using the maximum likelihood estimators. 
The range of the number of a/c desired on station at any given cycle will 
be set by the user. The number to be launched at any time is assumed 
equally likely within this range. This input parameter is a function of 
the estimated submarine density (i. e. , expected contact rate). The 
lower limit will be set at the number of a/c desired on station in the 
search (screening) phase, and the upper limit is set at the maximum 
practicable number of a/c to be launched during a multiple -contact 
phase. 

Briefly, the assumptions are: 

1. a/c will be relieved on station. 

2. Any desired length of operation can be set as an input. 

3. Major 120-hour checks will be completed prior to the operation. 

4. No resupply to the carrier is available. 

5. The launch-to-launch cycle for all ASW a/c is four hours. 

6. Minor maintenance, refueling, and rearming only can be 
performed on the flight deck. 
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7. Each maintenance spot is characterized with an independent 
exponential repair time with mean repair rate of X for 
around-the-clock operations. 

8. The number of a/c lost due to attrition is a Poisson random 
variable for each cycle period with parameter X (a/c 
accident/flying hours for a/c type). 

9. Any a/c lost by accident will not be returned to service due to 
either (a) physical loss at sea, or (b) insufficient mainte- 
nance capability aboard ship and lack of major parts. 

10. The number of a/c launched for each cycle is uniformly 

distributed between the upper and lower limits determined by 
the user. 



1 
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3. MODEL DESCRIPTION 



3. I The Theory 

A stochastic or random process is a collection of random variables 
indexed on some set T, (X (t) , t c T ) . In this case, time is the indexing 
set, and the Markovian assumption states that the future state of the 
process depends only on the state at the present time and not on its 
past history. Due to the cyclic nature of our problem, it is possible 
to increment time (T = (0, 1, . . . )) using the cycle time from launch 
to launch as the steps of unit time in a discrete Markov chain. It is 
assumed that the reader is familiar with the notion of a random variable 
as a function defined on a sample description space (S) on which the 
family of events or outcomes (E) of a probability function can be 
defined [4]. 

A discrete time Markov chain is described by a sequence of dis- 
crete valued random variables and is determined when the one -step 
transition probabilities of the state variables are specified, i. e. , a 

t 

conditional transition probability of a transition at time n for each pair 
of i, j = 0, 1, . . . , m (m being the number of states in the process) 
must be given. 

p..(n, n+ 1) = P [X(n + 1) = j | X(n) = i] 

If the transition probability functions depend only on the time difference, 
we have time homogeneity 
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p. . (n + 1 , 1) = p. . (0, 1) = p. . . 

1J 

The initial state of the system must be given either as a specific state 
or randomly as a probability distribution function over the possible 
states. 

The p.. (transition probabilities) are arranged in matrix form and 
satisfy: 

1. p. . £ 0 for 

y 

m 

2. Z p.. = 1 

j-o 1J 

for all i for 



~ 1 f ••• y , 

, i. e. , the rows of the transition matrix sum to 1 
the states within the description space [4]. 



20 



3. 2 The Model 



In order to establish the finite set of states (E) for the model, we 
shall consider two random variables defined as follows: 

(t) = The number of a/c flying at time t not having flown in 
the previous launch -to -launch interval. 

X^(t) = The number of a/c in or awaiting maintenance at timet. 

Now, we will consider the vector X(t) = [X^ (t) , X^ (t) ] as a pair of 
random variables and thereby have a bivariate stochastic process with 
the possible states ranging from (0, 0) to (A, N). 

0 ^ X^ (t) ^ A = No. of a/c desired on station, and 

0 ^ X^(t) ^ N = No. of a/c of given type aboard carrier. 

We will define an operating cycle as an interval unit of time. Process 
observations of X (t) will be made at successive unit interval launch 
times. To develop the p.. elements, consider a given time t for launching 
until A aircraft are flying or until the supply of ready a/ c is depleted. 
Those a/c failing the launch enter the maintenance state at this idealized 
point in time t (the total launching time required is much less than the 
total cycle time). At some time T, less than the launch-to-launch unit 
time interval, the a/c which were relieved on station return and land at 

the idealized point in time t + T. Some of these a/c will require mainte- 

nance and enter the maintenance queue. Those requiring only refueling 
and preflight inspection will enter a ready status to be tested for the next 
launch. 
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During the unit time interval, maintenance will be performed on 
those a/c in the not-ready status, and a certain number of aircraft will 
be repaired according to assumption 7. 

In summary, we start the system in some initial state (such as 

(0, 0) with no a/c flying or in maintenance) or start with a probability 

distribution Q(t ) over the states, E, at time t . We launch, recover, 
o o 

and repair a/c in the unit interval and repeat the process over each 
succeeding unit time interval until the end of the operating period. 

Knowing the transition probabilities within the unit time interval, we can 
develop the elements of the transition matrix, P, or {p 
These are the probabilities of going from the state of a a/c flying and 
i a/c in maintenance to (3 a/c flying and j a/ c in maintenance over the 
unit time interval. 

It was assumed in section 2 that A, the number of a/c to be launched, 
and N, the total number of a/c on board, are random variables, whereas 
they have been treated as constants so far in the development. To be 
analytically correct in including this feature, one should develop the 
appropriate quadr ivariate process. Such a development leads to too 
large a state space and the author chose to include these effects by using 
a Monte Carlo simulation technique. That is, at the beginning of each 
cycle, a random mechanism is used to determine the values on A and N. 

The probability of losing an a/c or changing the desired number to 
be launched is determined from the specified distributions at the beginning 



(p, j) } * 



'A 
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of each unit interval, and the resulting P matrix containing the 

p, .. is then recomputed. The probability distribution Q(t) 

(a, i) . (3 . j) 

over the states at any time t may be determined by the appropriate 
number of successive iterations of the Q vector times the P matrix, i. e. , 
Q(t) = P[X 2 (t) = p, X 2 (t) = j] = Q(t - l)xP . 

The probability of maintaining o' a/c on station over any given period 

of operation may be obtained at any unit time t (i. e. , the beginning of 
the next cycle) by summing out the appropriate maintenance state proba- 
bilities. Thus, P {<y a/c are flying at time t ) = 

Pr (X (t) = a ) = £ Pr (X (t) = a , X 2 (t) = i ) . 

i=o 

A mathematical comment appears to be in order. In the case of 
fixed A and N, the states of the Markov chain are positive recurrent; and 
steady-state probabilities can be found for the entire state space. In the 
case of decreasing N due to a/c attrition, this is not true; and (0, 0) 
becomes an absorbing state as time (t) goes to infinity. This latter 
consideration is not a realistic one for the operational period envisioned. 
Therefore, it is mathematically more feasible to use the former chain in 
conjunction with the Monte Carlo technique. 
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4. DEVELOPMENT OF THE TRANSITION MATRIX 



Perhaps the simplest way to view this development is to note the 
various transition probabilities incorporated in one -unit time cycle 
defined as follows: 

(1) y £ h = the launching transition probabilities at time t. This is the 

probability of taking f ready a/c, launching g successfully, and sending 
h into maintenance. Each a/c to be launched is considered a 
Bernoulli trial with probability of failure of p , which is estimable 
and subject to sensitivity analysis. The values of are: 

a. 0 if g > A , since only A a/c are desired; 

b. 0 if g 4- h > f ; it is impossible to launch and send into mainte- 
nance more a/c than are available; 

c. 0 if g < A, g + h < f ; launching continues until A a/c are flying 



or until all f are used up; 

f g f - or 

d. ( ) (1 - p ) (p ) if g < A, g + h = f , standard binomial 

& Y Y 

when all a/c in the ready state are used up but the A a/c are 

not launched; 
g + h - 1 g h 

e * ( Yi ) (1 - P) (p) if g = A, g + h> f , standard negative 
binomial for g successes in g 4* h - 1 trials. 



(2) n (m) = the landing transition probabilities which occur at time 
t 4- T . We must consider the probability that if there are a/c flying 
at time t then m a/c will enter maintenance at recovery time t4- T. 
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n a (m) will equal a standard binomial where p = the probability of 
equipment failure in flight: 



(3) p. . (t ) = the maintenance transition probabilities, i. e. , the proba- 
bility of repairing (i - j) a/c in time t . Two maintenance periods 
occur: the first starting at time t and ending at time t + T , the 
second starting at time t 4- T and ending at the end of the cycle, 

(t + 1) . Under assumption 7, the pulsed input, multiple exponential 
server queue is developed with D maintenance ’’spots" or servers 
each with identical, independent service rates, X . For each server, 

then, the probability of remaining occupied (given the server is 

” X. T 

busy) in time t = e . The probability of becoming free (i. e. , 

repairing an a/c) = 1 - e ^ T . The resulting queueing equations are: 



Three ranges of i (initial queue state) , j (final queue state) , and 
D become significant: 

a. When j ^ i ^ D, then not all spots are busy since there are 
fewer a/c in maintenance than spots. Each spot works inde- 
pendently; therefore, the solution to A is the binomial: 



rym) . 0(1 -p) 



ot - m 



(p) m , m = 0, 1 , . . . , a . 




B. dP. (t) / dt = -D\P. (t) + DXP. , , (t) 

i, n l, n l , n t 1 



for n ^ D . 



P..(t) = (j ) (1 - 



(i - j) 



- A.tj 



e 
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b. When D ^ j ^ i, then all spots are occupied throughout the total 
service time, and the closed form solution to B is the Poisson: 



p ij it) = 



(OX.)' 1 "-” e' DX * 

(i - j)i 



c. When j < D< i, then all spots are busy at the beginning of the 
service period, and some spots become idle during the service 
period. The explicit form solution of equation A is found using 
moment generating function transformation: 



D-l 



p ij <*> 



■ s. <?> (?) {(o^r) 



( i - D) - Un 
e 



- e 



- Xt 






k=o 




(The derivation of this solution is discussed in Appendix I. ) 



The figure on the following page will show the relationships of these 
transition probabilities within the unit time interval. 
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TRANSITION PROBABILITIES WITHIN THE UNIT CYCLE 

FIGURE 1 
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In order to develop each transition probability over the total unit 
time interval, we must consider all events taking place within the 
interval. Thus, to obtain the probability of going from cx a/c flying and 
i a/c in maintenance to (3 a/c flying and j a/c in maintenance, we start 
at the state (ex, i) at time t. At this time, a/c are launched and some 
1 a/c failing the launch enter maintenance. These i + 1 in maintenance 
are then serviced until time t + T when some k a/c are still in the 
maintenance state. At time t + T, of the a a/c previously flying, some 
m enter maintenance and (cx - m) enter the ready pool. Maintenance is 
continued on the (k + m) a/c for the remainder of the cycle (1 - T) , until 
the end of the unit time interval when j a/ c remain in the maintenance 
state. In functional form: 



N-cx-i i+1 cx 



P (or. i) > (P , j) 



r e s 



V N-a -i, p , 1 



l=o k=o m=o 
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5. SUMMARY 



Representative values for the mean repair rate and the landing and 
launching failure rates produced results in agreement with the sensitivity 
analysis by Collins on these parameters in [5], For failure proba- 
bilities less than . 5, and mean repair rate less than 12 hours, the effect 
of reducing the available maintenance time to 80% of the cycle time was 
negligible. Optimal loading and cycling policies can be determined for 
known values of these rates. 

The model affords the following checks: (1) the rows of each P 

matrix are summed as they are computed by the program; and (2) the 

probability distribution vector (Q J) is summed over the states. Each 

8 

summation was within 10 of one in the computer model. 

The user may substitute any available distribution over the interval 

of a/c desired on station. In order to keep A fixed, enter the desired 

value as both upper and lower limit ( A = ALOLIM = LUPLIM ) . For 

8 

fixed N, use a very small value for ALAM (such as 10 ). Subroutine 

KRAN is a uniform generator, using the half open interval (lower limit + 1, 
upper limit + 2) and a starting number as inputs. KRAN outputs an integer 
in this interval. Subroutine DRAW was used to provide some intuitive 
grasp of the results. DRAW was used in binary card form and is not 
essential to the main program. (The indicated associated statements 
must be removed, however. ) 
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The results of reasonable arbitrary parameter values, based on the 
author's experience, have shown that most of the probabilities concen- 
trate over a few states. Moreover ( , computation time increases rapidly 
as a function of N (no. of a/c), see Figure 2. This would indicate that a 
simple approximation to the model could be developed. One method 
presently being investigated to reduce computation time is to shrink the 
probability state space to include only tholse significant states and, thus, 
reduce the size of the transition matrix. Alternatively, the eigenvector, 
eigenvalue representation of the P matrix, might be used. 

Originally, it was hoped to utilize the data from the Fleet ASW Data 
Analysis Program (FADAP) to attempt a verification of the model with 
its real-world counterpart. The only method available at present for 
obtaining the necessary data is by direct observation or a program of 
data collection, as suggested by Collins [ 5] . 

Many fruitful areas of investigation exist: 

(1) Attrition has been simply modeled by the Poisson method. The two 
components of attrition, accidents and supply shortage, can be more 
accurately modeled and used to develop logistic schedules for 
maintenance and supply. One simple technique is to assume each 
component is independent and Poisson, and estimate a supply failure 
rate for AOCP attrition from past data. With these assumptions, 
the total attrition is Poisson, with the parameter equal to the sum 
of the accident and supply failure rates. 
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(2) The model can be modified to make the number of maintenance spots 
available for any cycle a variable function of time, D(t) * 

(3) An investigation of the Markovian assumption validity as the cycle 
times become smaller and smaller. 

(4) Development of a continuous time model. 

(5) Modification of the model to simulate resupply by COD. 

(6) A study of the distribution of submarine contacts to determine the 
validity of the uniform a/ c demand assumption. 
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PROGRAM ASSEMBLY AND COMPUTATION TIME 
FOR ONE TRANSITION MATRIX (P) AS A FUNCTION 
OF THE TOTAL NUMBER OF AIRCRAFT (N) 



FIGURE 2 
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APPENDIX I 



EXPLICIT SOLUTIONS OF THE 
MAINTENANCE QUEUEING EQUATIONS 



The queueing equations for the pulsed input queue are essentially 
the pure death process given in [ 1 ] and [4] as problems and developed 
by Collins in [5], The equations are: 

dP (t) 

A. — = - n X P. (t) + (n + 1) X. P. ,(t) forO<n<D 

dt l, n i, n + 1 

dP i n (t) 

B. = - nXP. (t) + DXP. (t) for n ^ D 

dt i, n i, n + 1 



where P. . (0) = A. . and P. . (t) = 0 for i < j , since no input (arrivals) 
occur during the service time. 

Equation B is solved directly in closed form: 



P i n (t) = 

i, n 



(D X t) (n - l) e- XDt 
(n - i) ! 



Now transforming the first equation (A) using the moment gener- 
ating function (MGF), r 

' ^ 

D - 1 _ 

G(s, t ) = L s n p (t) , 

^ n 

n=0 

as outlined in [4] (Chapter 7), and its partial derivatives: 
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( 1 ) 



dG 



D-l 

dt £ S P n (t) 

dt n=0 



( 2 ) 



dG ^ n “ 1 D / f \ 

— ; — = Zj ns P (t) 

ds n 

n=0 



Where P (t) denotes the conditional probability P. (t) , by substituting 
(A) into (1), properly identifying the first summation with (2), and 
changing the second summation index to r = n + 1 , we get: 



dG 

dt 



dG 



D 



X s — — + X L r s 

ds 

r = 0 



r - 1 



P (t) , or 
r 



< 3 > IF = 'Ms-1) ig- + XDs D - 1 P D (t) , 



since 



£ r s r " 1 P (t) = + Ds D ’ 1 P (t) 

^ r as D 

r =0 



Next, replace the partial differential equation (3) with a system of 

>y 

ordinary differential equations using the Lagrangian auxiliary equations: 



ds 



dt 



1 X (s - 1) 



dz 



X D s D " 1 P D (t) 



The solution to the first equation (using the first two differentials) 



is: 



Xt = ln(s-l) + C' 
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and hence 



^ ^ t 

s = e +1 
or 

_ , — \ t . . . 

= e (s-1) 

The second equation is: (using first and third differentials) 

dz = - \D (C e Xt + 1) D ' 1 P (t) dt . 



Using the solution to (B) where m = i - D to replace P (t) and 



integrating, term wise, the binomial expansion of ( e ^ + 1 ) ^ 



z = 



UP) 

m! 



m + 1 D _i 



jfo (VJ C ‘ j J* 



m e -\( D -j)t dt 



where the integral is evaluated as: 



m 

- X 



t k e -MD-j)t 



k=0 ( X (D - j) ) 



m - k + 1 



m! 

“kT + °2 



Thus , 



C 2 = z + e 



, _ D-l m k 

-XDt ^ , D i „ ( \ D t) 



s (?) (s - l) J X 

j=0 k=0 



k! 



(A) 



m - k 



and the general solution is 0 (C^ , ) , where 0 is an arbitrary function 

and 

C 1 = u (s, t, z) 

and 



C 2 = v (s, t, z) . 
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To get our particular solution, use the boundary conditions for G(s, t) 



( 1 ) for s = 1 , 



D-l 

G ( 1 , t) = 2 P (t) 

^ n 



= Pr [no. in maintenance at t is < D | i at t = 0 ] 



G ( 1 , t) 



i-D 

2 

n=0 



.- tD1 |imr 

n ! 



- t^t) 



where 



u ( 1 , t, z) = = 0 



v (1 , t, z) = C 2 



- X Dt 

z + e 



m 

2 

k=0 



( X Dt) k 
k! 



so 



C 2 = z + ^ i W 



(2) for t = 0 , 

G ( s, 0) 



where 



D-l 

2 

n=0 



s n p (0) 

n 



0 , since i ^ n > D 



u (s, 0, z) = (s - 1) 



V (s, 0, z) 



D-l 
z + 2 

j=0 



(?) 



• ‘> j CdTt) 



m 



Thus, 



G (s, 0) 



z 



D-l 
+ 2 
j=0 



(?) 



J G^t) 



m 



C 



2 
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Substituting the general value for C 



2 



above : 



D-i n . . . 

G (s, t) = 0 (u, v) = L ( : ) (s - 1) J e ' K J 

j = 0 



CdTj) 



m 



D £ ^ ,-XD. 

j=0 j k=0 k ' 



(^l) 



m 



Rearranging terms, 



D-l D- 1 



g(s, t) = s s n z (j> (° ) (-n j r 

n=0 J=n " 1 L V. D - 



m 



Xt j 



- X Dt ™ ( X Dt) k 

© 2^ 



k=0 



k! 



(dTt) 



m - k 



where P (t) = the coefficient of s 11 . 

n 
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APPENDIX II 



THE LOGICAL F' OVf DIAGRAM OF THE COMPUTER PROGRAM 

( A,F,p , PGAM TIME, ICYCLE, JCVCLE 
7 LAM, ALAM, AUPLIM, ALOLIM, ETC.) 



\ SET ; 

INPUTS 




COMPUTE ( ra ) =PRFUA ( I AA , MI ) 



4 

v? 

J 

COMPUTE Pij(t) = PTRM(I,IJ,IT)| 
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APPENDIX III 



THE COMPUTER PROGRAM 
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nnnnnnnnnnnnnononoon 






-B?NARY l ;56? AN,0M9/S/1S/2S/E/45 = 54 ’ 15 » 30000 ’ 5 * 

( RELOCOM. 

-FTN »E. 

PROGRAM MARKOV 

THIS PROGRAM IS A NONSTATIONARY BIVARIATE 
OPERATIONS. THE RANDOM VARIABLES ARE THE 
BEGINNING OF ANY GIVEN LAUNCH CYCLE. THE 
THE MODEL IS 16(NA). THE RANGE OF A/C TO 

INTERVAL IS 0 TO 6 A/C. THE FOLLOWING INPUTS ARE REQUIRED. 
ID- THE NO. OF INDEPENDENT MAINTENANCE SPOTS 
NA= TOTAL NO. OF A/C TYPE ON BOARD 



809 



T IME = T IME 
FLAM=MEAN 



FROM LAUNCH 
REPAIR TIME 



TO RFCOVERY /LAUNCH 
PER SPOT/LAUNCH TO 



000 
000 
000 
000 
000 

MARKOV CHAIN MODEL OF ASW A/C 000 
NUMBER OF A/C FLYING AT THE 000 
MAXIMUM NO. OF A/C ALLOWED IN 000 
BE LAUNCHED AT ANY GIVEN 000 

000 
001 ' 
001 



TO LAUNCH CYCLE T I ME ( HRS ) 00 1, 
LAUNCH CYCLE TIME (HRS)OOl. 



n GA ^o« PR0BAR 1 L 1 TY FILING LAUNCHIM.L.EST. FROM PAST DATA) 

PROBABILITY OF A/C FAILURE DURING FLIGHT REQUIRING MAINTENANCE 
AT LANDING (M.L. ESTIMATOR FROM PAST DATA) 

QI- THE PROBABILITY DISTRIBUTION VECTOR OVER ALL POSSIBLE STATES 
(7 X 17 = H9) SUCH THAT THE SUM OF ALL QI ( I ) = 1. j H IS 

ICVfLF AND INPUTTED By USING A DATA STATEMENTOO 1' 

ICYCLF - NO. CYCLES DESIRED FOR OPERATION nn? , 

= LAUNCH TO LAUNCH T I ME ( HRS ) ( TOT . T IME= I CYCLE X JCYCLE) 

ACCIDENT RATE FOR TYPE A/C (ACCIDENT/HOURS)' JLYCLt > 

- DESIRED LOWER LIMIT ON A 
= DESIRED UPPER LIMIT ON A 



001 - 
001 ' 
00 1 < 
oor 
0011 



JCYCLE 
ALAM = 

ALOL IM 
AUPL IM 

COMMON FLAM.TIMF 
TYPE DOUBLE' FLAM 
COMMON PTRM, GAMMA, PR ,PRFMA , ID 
DIMENSION BC( 17) ,A( 17) ,FBC ( 17) 

DIMENSION PTRM(17»17»2) , GAMMA ( 17,7,17) »PRFMA(7»7) 
DIMENSION PR( 119,7,17) ,QI (119) ,QJ( 119) 

DIMENSION FJPLOT ( 119) ,JT(12) 

ENTER DATA CARDS HERE 

DATA ( (QIC I) ♦1 = 1,119) = . 2, 1 6 ( . 05) *102( .0) ) 

NA= 16 
ALAM= .01 
ID = 8 
FLAM=3.0 
PGAM=P=.4 
IYY = 13421773 
TIME = .125 
ICYCLE=20 
JCYCLE=4 
ALOL IM=4. 

AUPL I M = 6 • 

END OF DATA CARDS 
AL=ALOLIM+l. $ AU=AUPLIM +2. 

UN I TT = 1 • 

N=NA+1 
I AMAX= 7 
I A L A S T = 0 
D=FLOATF(ID) 

NLASf=NEWN=N 
KT= 1 

I A = KRAN ( AL , AU , IYY) 

IF(KT.EQ.l) 113,115 



002! 
00 21 
0021 
002; 
002: 
00 2 * 
002f 
002C 
0021 
00 26 
0025 
003C 



0031 

0032 

0033 

0034 

0035 

0036 

0037 
00 3§ 

0039 

0040 

0041 
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115 T1=-L0GF( .000000001 + RANF ( -1 ) ) *2 . 30258/ALAM 0042 

IF(Tl-TFLC) 130,131 *132 0043 

130 T2 = -LOGF( .000000001 + RANF ( -1 )) *2.. 30258/ALAM 0044 

I F ( T1+T2-TFLC) 230, 231, 1M 0045 

230 T3=-LOGF( .000000001 + RANF ( -1 ) ) *2 . 302 58 / ALAM 0046 

I F ( T1+T2+T3-TFLC) 331,331,231 0047 

331 NEWN=NLAST-3 $ GO T0113 0048 

231 NEWN=NLAST-2 $ GO T0113 0049 

132 NEWN=NLAST $ GO T0113 0050 

131 NEWN=NLAST- 1 0051 

113 PRINT 8882 ♦ I A ,NEWN 0052 

IF(NEWN-IA) 15,13,13 . 0053 

15 I A=NEWN 0054 

13 IF(IALAST) 11,12,11 • 0055 

12 CONTINUE 0056 

FROM THIS NEXT STATEMENT TO NO. 483 IS CONCERNED ONLY WITH THE GRAPH 0057 
DO 482 1=1,12 0058 

482 JT ( I )=8H 0059 

JT ( 1 )=8HE( A/C) = 0060 

JT ( 3 )=8HSPOTS = 0061 

JT ( 5 ) = 8H T = 0062 

JT ( 7 ) =8H J VS QJ 0063 

JT ( 8 )=8HVECT0R 0064 

JT ( 9 ) =8H N = 0065 

JT ( 1 1 ) =8H A = 0066 

DO 483 1=1,119 ^ 0067 

F I = I 0068 

483 FJPLOTC I )=FI 0069 

I ALAST = I A 0070 

DO 1235 1=1,17 0071 

DO 1235 J= 1 , I AMAX * 0072 

DO 1235 K= 1 » 17 0073 

1235 GAMMA ( I ,J,K)=0.0 0074 

AT THIS PT THE LANDING TRANSITION PROBABILITIES ARE COMPUTED 0075 

DO 300 I AA= 1,1 AMAX 0076 

DO 301 M I = 1 ♦ I AMAX 0077 

IF{ IAA-MD31, 32*33 0078 

31 PRFMAf IAA,MI)=0. 0079 

. GO TO 301 0080 

32 MM1=MI— 1 0081 

PRFMA( IAA,MI ) =P**MM1 0082 

GO TO 301 0083 

33 I AMI = I AA- 1 0084 

‘ MM1 =M I - 1 0085 

BC ( 1 ) = 1 .0 0086 

PROD=FLOATF( IAA-MI ) 0087 

DO 50 I P = 2 ♦ M I 0088 

AIP=FLOATF( IP-1 ) 0089 

PROD = PROD +1.0 0090 

50 BC( IP)=PROD*BC( IP-1) /AIP 0091 

I GO= I AA-M I 0092 

PRFMA ( I AA ,M I ) = ( BC ( M I ) * ( 1 . 0-P ) ** ( I GO ) ) *P**MM 1 0093 

361 eOMTfNME 

300 CONTINUE 0095 

AT THIS PT THE MAINTENANCE TRANSITION PROBABILITIES ARE COMPUTED 0096 

DO 100 I T = 1 ,2 0097 
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IF( IT-1 )25,25 ,26 

25 TAU = TIME 
GO TO 28 

26 TAU = UNITT-TIME 
28 DO 101 1=1, N 

DO 102 I J = 1 , N 
IF (I-IJ) 14,199,17 
199 IF(I-ID) 19,19,1999 
1999 PTRM( I , IJ»IT)=EXPFC-FLAM*TAU*D) 

GO TO 102 

14 PTRM(I,IJ,IT)=0. 

GO TO 102 

19 FJM1=FL0ATF ( I J-l ) 

PTRM( I ,IJ,IT)=EXPF(-FLAM*TAU*FJM1) 

GO TO 102 

17 IF ( I— ID— 1 ) 1,1,2 

1 BC ( 1 ) = 1 .0 
PROD=FLOATF( I-IJ) 

DO 10 IP =2 , I J 

A I P = FLO ATF C I P — 1 ) 

PROD = PROD + 1.0 
10 BC ( IP) =PROD*BC( IP-1J/AIP 
ELT=EXPF(-FLAM*TAU) 

PTRM ( I »IJ, I T ) = BC ( IJ) * ( 1 ,-ELT )**( I-IJ)#ELT**( IJ-1 ) 

GO TO 102 1 1 

2 IF(IJ-l-iD) 22,24,24 
22 CONTINUE 

CALL P IDU , I J , IT ) 

GO TO 102 
24 D=FLOATF ( ID ) 

ELDT = FXPF ( -D-x-FL AM*T AU ) 

FACT = 1.0 
A( 1 ) =1 .0 
MM = I-IJ 
DO 20 M=2 , MM 
FACT = FACT + 1.0 
20 Af M ) =A( M-l ) *FACT 

201 PTRMf I »IJ»IT) = (D*FLAM*TAU)**(I-IJ)*E|_DT/A( I - I J i 

102 CONTINUE . /A ‘ 1 1J) 

101 CONTINUE 
100 CONTINUE 
GO TO 120 
11 CONTINUE 

IF( IA-IALAST) 120,121,120 
121 IF(NEWN-NLAST) 111,117,111 
117 I ALPH= I ALASTM1 $ GO TO 801 
120 CONTINUE 

AT 00 I 2<H°'fF=UN LAUNCH,NG TRAN SIT10N. PROBABILITIES ARE COMPUTEO 

IGM = XMINOF ( IA, IFF) 

DO 203 IG= 1 , I GM 
I GM 1=1 G-l 
G>0 202 IH = liN 
I HM1 = I H-l 

BPROD= ( ( l.-PGAM )**IGM1 ) *( PGAM** IHM1 ) 

86 IF(IG-IA) 91,87,84 



OC 
OC 
01 
or; 
01 
oi ! 
01 
or 
01 
oi: 
or 
oi 
oi 
oi 
oi. 
01 ' 
o i :j 
oi : 
oi : 
oi: 
on 

Oil" 

oi; 
012 
012 
012 
012 
012 
012 
012 
012 
012 
013 
013 
013 
013 
013' 
0 1 3 : 
0 1 3c 
0 1 3' 
01 3fc : 
013' 

0 1 4( 
014! , 

014; : 

014; 

0144 

0 1 4f 

014« 

0141 

014? 

0149 

015g 

6151 

0152 

0153 
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91 IF< IG+IHM1-IFF) 84,82,84 0154 

87 IF( IG+IHM1-IFF)85, 85,84 0155 

84 GAMMA ( I FF » I G » IH ) =0 • 0156 

GO TO 202 0157 

82 BC< 1 )=1.0 0158 

PROD=FLOATF( IFF-IG) 0159 

DO 30 I P=2 » IG 0160 

AIP=FLOATF( IP-1 ) 0161 

PROD = PROD + 1.0 0162 

30 BC ( IP) =PROD * BC ( I P- 1 ) /A I P 0163 

I HM1= I H-l 0164 

TEMP= PGAM#*IHM1 • 0165 

TEMP1=( l.-PGAM)**IGMl 0166 

BPROD = TEMP*TEMP1 0167 

GAMMA { IFF,IG»IH)=BC( IG)*BPROD 0168 

GO TO 202 0169 

85 FBC ( 1 ) = 1 • 0 0170 

PROD=FLOATF ( IGM1-1 ) 0171 

DO 40 I P=2 ♦ I H 0172 

A I P = FLOATF ( IP-1 ) 0173 

PROD = PROD +1.0 0174 

40 FBC ( IP) =P ROD* FBC ( I P- 1 ) / A I P 0175 

GAMMA ( I FF , I G » IH)=FBC( IH)*BPROD 0176 

202 CONTINUE 0177 

203 CONTINUE 0178 

204 CONTINUE ' 0179 

REMOVE CARDS FROM HERE TO NO 999 IF PRINT OUT NOT DESIRED 0180 

PRINT 9*( ( ( I*IJ»IT»PTRM( I»IJ»IT)»IT=1»2) »IJ=1»N)»I=1»N) 0181 

9 FORMAT ( 1H 1 / ( 2 ( 6H PTRM ( I 2 » 1H , I 2 ♦ 1H , I 2 , 3H ) = E 14. 5 ) ) ) 0182 

PRINT 99 » ( ( ( IFF, IG»IH,GAMMA( IFF»IG, IH) ♦ I FF = 1 , N ) » I G= 1 » I A ) »IH=1,N) 0183 

99 FORMAT ( 1H1/ ( 2 ( 7H GAMMA ( I 2 » IH , I 2 » IH , I 2 , 3H ) = E14.5))) 0184 

PRINT 999, ( ( IAA,MI ,PRFMA( I AA,MI ) , I AA = 1 , IAMAX ) ,MI = 1 , I AMAX) 0185 

999 FORMAT ( IH 1 / ( 2 ( 7H PRFMA ( I 2 » IH , I 2 ,3H ) = E 1 4. 5 ) ) ) 0186 

NOW THE TRANSITION MATRIX MUST BE ZEROED 0187 

111 CONTINUE i 0188 

DO 899 J= 1 » 1 19 0189 

DO 899 K = 1 ,7 0190 

DO 899 L= 1 » 17 0191 

899 PR( J,K,L) =0.0 0192 

START COMPUTING THE ELEMENTS OF EACH ROW, I=IX+ (ALPHA - 1 ) X TOTAL- A/C 0193 
DO 1000 I ALPH = 1 , IALAST 0194 

801 CONTINUE 0195 

DO 1100 I X= 1 ♦ NLAST 0196 

COMPUTE THE P ELEMENTS OF THE IAPH,IX ROW AND SUM THE ROW 0197 

T $UM=0 . 0198 

I=IX+( IALPH-1 )*N 0199 

DO 800 I BET A= 1 ♦ I A .* ' 0200 

RSUM=0 . 0 0201 

DO 900 I Y = 1 ,NEWN 0202 

PR ( I , I BET A, I Y ) = 0. 0203 

I L I M=NEWN- I ALPH- I X+2 0204 

PSUM=0.0 0205 

SUM-0.0 0206 

SUML=0.0 0207 

DO 500 I L= 1 ♦ I L I M 0208 

KL I M= I X+ I L- 1 0209 
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1 



701 

702 



700 

600 

500 

900 

800 

888 

1100 

1000 



IXPIL = IX + I L - 1 
SUMM=0 • 

DO 600 M I = 1 , I ALPH 
SUMK=0 • 

DO 700 IK=1»KLIM 
IKPMI = IK +MI-1 
I F ( IXPIL-NEWN) 701,701,700 
IF( IKPMI-NEWN) 702,702,700 
GAMH=GAMMA ( I L IM , I BET A , I L ) 
PTRMH1 = PTRM( IXPIL, IK, 1 ) 
PRFMAH = PRFMA ( I ALPH ,M I ) 
PTRMH2 = PTRM( IKPMI, IV, 2) 

SUM = GAMH * PTRMH1 * PRFMAH 
SUMK =SUMK + SUM 
PSUM=PSUM+SUM 
CONTINUE 

SUMM = SUMM + SUMK 
CONTINUE 

SUML = SUML +SUMM 

PSUM2 =SUML 

CONTINUE 

RSUM=RSUM+PSUM 

PR ( I , I BETA, I Y ) =PSUM 

CONTINUE 

TSUM=TSUM+RSUM 

CONTINUE 



* PTRMH2 



TSUM, IALPH, IX 
7H 



TSUM = E15. 8*215) 



HERE TO 889 IF P MATRIX PRINT OUT NOT DESIRED 



889 

890 

898 
C NOW 
805 
807 



CAT TH 



PRINT 888 , 

FORMAT ( 

CONTINUE 
CONTINUE 

REMOVE CARD FROM 
DO 889 J-1,17 
DO 889 K = 1 » 7 
D0889 L = 1 , 17 
I = J+ ( K- 1 )*N 

PRINT 890 , ( PR ( I ,LP ,L ) ,LP = 1 , IAMAX) ,K, J,L 
FORMAT(7E14.5,2HJ=I2,5HK=l,A,2HL=I2) 

DO 898 1=1,119 
Q J ( I ) =0 .0 

MULTIPLY Q I AND P TO GET QJ 
PRINT 807,KT, IALAST, IA 

FORMAT ( 1H1,13HQ VECTOR CASE 13/// 15,15) 
DO 802 IBETA = 1 ,7 
DO 902 I Y = 1 , 1 7 
IS POINT CALCULATE THE 
J= I Y + ( I BETA-1 ) *N 
0P1=0. 



(IBETA,IY)TH ELEMENT OF THE QJ VECTOR 



OP = 0. 

DO 2001 I ALPH= 1,7 
DO 2201 I X= 1 » 17 
I = IX+('I ALPH-1 )*N 
QP 1 =Q I ( I ) *PR ( I » I BETA » I Y ) 
r5P = QP-f©Pi 
2201 CONTINUE 
2001 CONTINUE 
Q J ( J ) =QP 



021 
021 
021 'J 
021 
021 
021 
021; 

02 l: : 

0 21, 

0 2 1 f ;t 
02 2( , 
022: 
022; 
022: 
0224 ' 
02 2f 
022e 
02 21 ■ 
02 2i 
02 2S 
023C 

0231 

0232 

0233 
02 34 

0235 

0236 

0237 

0238 
0239,. 

0240 

0241 

0242 

0243 

0244 

0245 

0246 

0247 

0248 

0249 
0 250 

0251 

0252 
02 53 

0254 

0255 

0256 

0257 

0258 

0259 

0260 
0261 
0262 

0263 

0264 

0265 
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PRINT 8882 , I BETA, IY, J,QP 0266 

•882 FORMAT! 214, 4H QJ(I3,3H )= E14.8) 0267 

90 ? CONTINUE 0268 

802 CONTINUE 0269 

[CHECK THE SUM OF THE Q VECTOR 0270 

QSUM=0 • 0271 

DO 808 J= 1 ♦ 1 1 9 0272 

808 0SUM=0 J ( J ) +QSUM 0273 

PRINT 8883, OSUM 0274 

3883 FORMAT ( 6H QSUM= E15.9) 0275 

DO 333 I = 18,119 0276 

K = (I-D/17 • 0277 

FK=FLOATF(K) 0278 

FMEAN= FK*QJ ( I ) +FMEAN 0279 

333 CONTINUE 0280 

TFLC=FMEAN*FLOATF( JCYCLE) 0281 

PRINT 335 , FMEAN 0282 

335 FORMAT! 17HMEAN A/C FLYING = F10.4) 0283 

STATEMENTS FROM THIS POINT TO THE CALL DRAW STATEMENT REFER TO GRAPH 0284 
JT ( 2 ) = I CODE (FMEAN ) 0285 

JT ( 4 ) = I CODE ( D ) 0286 

FKT=FLOATF(KT) 0287 

JT ( 6 ) = I CODE ( FKT ) 0288 

FN=FLOATF(NEWN-l) 0289 

JT! 10) = ICODE! FN ) 0290 

FIAA=FLOATF( IA-1 ) 0291 

JT ( 12)=IC0DE(FIAA) . 0292 

CALL DRAW! 119, FJPLOT,QJ,0,0,4H , JT ,0 , 0 , 0 , 0 , 0 , 0 , 8 » 8 , 0 , L AST ) 0293 

FMEAN = 0. 0294 

NEXT WE MUST MULTIPLY OJ AND P TO G p T QK AND SO ON • • • ( QK+. . . N ) 0295 

KT=KT+1 0296 

IF(KT-ICYCLE) 803,803,806 0297 

803 DO 804 1=1,119 0298 

804 OKI ) =QJ ( I ) 0299 

IALASTM1=IALAST 0300 

I ALAST = I A 0301 

NLAST=NEWN 0302 

GO TO 809 0303 

806 STOP 06 0304 

END v 0305 

SURROUTINF PID(T,J,IT) 0306 

COMMON FLAM, TIME 0307 

COMMON PTRM, GAMMA ,PR ,PRFMA , I D 0308 

TYPE DOUBLE BC,BDC,PROD , D I D3 , D I D4 , D I D5 , DEXP 0309 

TYPE DOUBLE DAN , D I D 1 , D I D2 , SUM , DN , ANM 1 , FAC , COF , PSUM , PTR , FLAM , TAU , D 0310 
DIMENSION PTRM ( 1 7 » 1'7 » 2 ) , BC (11) ,BDC(11) 0311 

DIMENSION GAMMA! 17,7,17) , PRFMA(7,7), PR(119,7,17) 0312 

D=FLOATF(ID) 0313 

I DP 1 = I D + l 0314 

IF! IT-1 )25, 25, 26 0315 

25 TAU = TIME 0316 

GO TO 28 0317 

26 TAU= 1 .-TIME 0318 

28 CONTINUF 0319 

I MDP1= I - I D 0320 

PTR=0.0 0321 
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nnnn.nnnnnnnnn 



PSUM=0. 032 

DO 200 NJ = J, ID 032 

C DEVELOP N TAKEN J AT A TIME AND D TAKEN N AT A TIME 032 

BC ( 1 ) = 1 . 0 032 

PROD=FLOATF(NJ-J) 032 

DO 10 I P = 2 * J 032 

AIP=FLOATF( IP-1 ) 032 

PROD = PROD+l.O 032 

10 BC ( I P ) =PROD* BC ( I P-1 ) / A I P 033 

BDC ( 1 ) = 1 • 0 ' 033 

PROD=FLOATF ( IDP1-NJ) 033 

DO 20 I Q=2 * NJ ' 033 

AIQ = FLOATF( IQ-1 ) 033* 

PROD=PROD+1.0 033' 

20 BDC( IQ) =PROD*BDC( IQ-1 ) /AIQ 033 

C0F=BC( J)*BDC(NJ)*(-1)**(NJ-J) 033 

ANM1=FL0ATF( NJ-1) 033 

DAN=D/ ( D-ANM1 ) 033 

DID4=DFXP(-FLAM*T AU* AN Ml ) 034 

DID1=(DAN**( I-IDP1) )*DID4 034 

SUM=0. 034 

DN=0. 034 

DO 201 K= 1 * IMDP1 034 

FAC= 1 • 034 

KM1 =K-1 • 034 

PROD=0 • - ‘ 034 

DO 11 I K = 1 * KM1 034 

PROD = PROD+ 1 • 034 

11 FAC=FAC*PROD 03! 

I M I DK= I - 1 D-K 031 

SUM=( ( FLAM*D*TAU)**KM1 )*DAN**IMIDK / FAC 03! 

201 DN = DN+SUM 03! 

DID3=DEXP( -FLAM*D*TAU ) 03! 

DID2=DN*DID3 03! 

DID5 = DID1-DID2 03! 

PSUM=COF*DID5 03! 

200 PTR =PTR +PSUM 03! 

103 CONTINUE 03! 

PTRM( I * J» IT ) =PTR 03i 

102 CONTINUE 03i 

101 CONTINUE 03i d 

END 03' 

FUNCTION KRAN ( A » B * I Y ) 03 

03 

THIS ROUTINE RETURNS AN UNIFORMLY DISTRIBUTED RANDOM INTEGER 03 

03| 

THIS ROUTINE RETURNS A INTEGER RANDOM NUMBER .GE. TO A 03 

.LT. B 03 

A = BOTTOM LIMIT (INCLUDED) FOR THE RANDOM NUMBER 03 

B = TOP LIMIT (NOT INCLUDED) FOR THE RANDOM NUMBER 03 

SET IY ONLY ONCE IN MAIN PROGRAM FOR EACH SET OF RANDOM NUMBERS 03 

SOME GOOD STARTING VALUES FOR IY FOLLOW 03 

13421773 03 

33554433 03 

8426219 03 

42758321 03 
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56237485 

62104023 

ANY OF THESE MAY BE USED 

THIS ROUTINE MAY BE USED IN FORTRAN 60 OR 63 
IY = 3125 * IY 

IY = IY - ( IY/67108864) * 67108864 
FY = IY 

- KRAN = FY/67108864. * (B-A) + A 
RETURN 
END 

FINIS 
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APPENDIX IV 



SAMPLE RESULTS 

The following pages present the values of the elements of the proba- 
bility distribution vector (Q J) and its graphical plot for five consecutive 
iterations, i. e. , Q x P n for n = 1, 2, . . . , 5 . The inputs are those 
shown on the first page of Appendix III between statement No. 30 and 
No. 31. The printouts of the transition matrices and their computational 
elements are omitted. The plot was made using the DRAW subroutine in 
the U. S. Naval Postgraduate School computer facility library. Each 
vector printout contains the values of all 119 states possible (7 x 17) and 
is headed by the past value of A + 1 and the next value of A + 1. The two 
indices preceding each element represent 3+1 and j + 1, in the notation 
of section 3. For example, in the first row on the next page, the n l 1" 
indicates that the probability of being in state (0, 0) after one iteration 
is = . 026, where the value of A is 4 over the first iteration. Each graph 
is labeled with the expected value of a/c flying, the number of mainte- 
nance spots available, the vector number (T), total number of a/c avail- 
able (N), and the desired number of a/c on station (A). The M E n notation 
indicates the power of 10 to multiply by. This sample run demonstrates 
the loss in total a/c and variable a/c on station. 
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